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ABSTRACT 

Late in the gaseous phase of a protostellar disk, centimeter-sized bodies probably 
settle into a thin "dust layer" at the midplane. A velocity difference between the dust 
layer and the gas gives rise to turbulence, which prevents further settling and direct 
gravitational instability of the layer. The associated drag on the surface of the layer 
causes orbital decay in a few thousand years — as opposed to a few hundred years for 
an isolated meter-sized body. Within this widely-accepted theoretical framework, we 
show that the turbulent drag causes radial instabilities even if the selfgravity of the 
layer is negligible. We formulate axisymmetric, height-integrated dynamical equations 
for the layer that incorporate turbulent diffusion of mass and momentum in radius and 
height, vertical settling, selfgravity, and resistance to compression due to gas entrained 
within the dust layer. In steady-state, the equations describe the inward radial drift of 
a uniform dust layer. In perturbation, overdense rings form on an orbital timescale with 
widths comparable to the dust-layer thickness. Selfgravity is almost irrelevant to the 
linear growth rate but will eventually fragment and collapse the rings into planetesimals 
larger than a kilometer. We estimate that the drag instability is most efficient at 1 AU 
when most of the "dust" mass lies in the size range 0.1-10 meters. 

Subject headings: planetesimals; planetary formation 



1. INTRODUCTION 

Meteoritic evidence and theoretical considerations indicate that planet formation begins with 
collisional agglomeration of dust particles into larger bodies during the lifetime of a gaseous proto- 
stellar disk (< 10 7 yr). 

Because the gas is partly supported against the gravity of the central star by pressure, the 
orbital velocity of the disk is slightly lower than the keplerian circular velocity (Vk = \J GM/r). If 
c s is the sound speed of the gas, then 

_ Vk - V gas ^ r_8P (c_ 

71 ~ V K ~ 2p g Vj2 dr ~ V^k 
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Solid bodies bodies couple to the gas by aerodynamic drag and orbit at intermediate velocities. 
Small grains move almost exactly with the gas, while large planetesimals move nearly at Vk- In 
both limits, drag causes very little dissipation of orbital energy. Dissipation is maximal for bodies 
of intermediate size. At 1 AU in a minimum-mass solar nebula (Hayashi et al. 1985), the orbital 
decay time is only a few hundred years for particle radii r p ~ 1 m (cf. §2). Therefore, the processes 
responsible for growth must move bodies rapidly through the intermediate size range. Since the 
orbital decay time scales roughly linearly with r p above ~ 1 m, planetesimals larger than r p ~ 10 km 
experience little decay during the lifetime of the nebula. 

While collisional agglomeration may be sufficient to reach > km size before the orbits decay, 
collective instability is an attractive alternative. Safronov (1969) and Goldreich Sz Ward (1973, 
henceforth GW) suggested that planetesimals may form by gravitational instability. As an initially 
turbulent disk becomes quiescent, or as the solid bodies grow to radii r p ~ 1 cm, the particles settle 
into an increasingly dense layer — the "dust layer" — at the midplane. When the volume-averaged 
density of the dust layer approaches p* = 3M/ Aira 3 , where M is the mass of the star and a is the 
orbital radius, the layer becomes gravitationally unstable. Self-gravity then fragments the layer 
into planetesimals with radii ~ 1 km. 

Gravitational instability along the above lines is now regarded as unlikely. The dust layer 
would have to be very much thinner than the gas disk, since the mass fraction of solids at solar 
abundances is small (~ 10~ 2 ), and since the gas density itself is already ~ W~ 2 p* in a minimum- 
mass solar nebula. On the other hand, the velocity shear between the dust layer and the gas 
probably causes turbulence, even if there is no other source of turbulence in the gas disk. Although 
GW called attention to this "boundary-layer" turbulence, they considered its implications only for 
the exchange of momentum between the dust layer and the gas. Weidenschilling (1980) pointed 
out that turbulence will also tend to limit the density of the dust layer by mixing. More recently, 
Cuzzi et al. (1993, henceforth CDC) have confirmed this with a detailed analysis particle settling, 
shear-driven turbulence, and turbulent vertical mixing. These authors conclude that gravitational 
instability is unlikely to occur at 1 AU until r p 3> 1 m, if at all. 

Even if direct gravitational instability does not occur, there may be secular instabilities — that 
is, instabilities relying upon dissipation. GW mentioned the possibility of secular effects but did 
not elaborate. They may have had in mind the following axisymmetric instability: Consider a 
steady background state in which the surface mass density of the dust layer (S p ) is uniform on 
scales <C r and drifts radially inward at a rate proportional to the drag against the gas. Let there 
be a positive perturbation in S p on a ring of radial width (Ar) smaller than the thickness of the 
gas layer (H g ) but larger than that of the dust layer (H p ). The (slight) selfgravity of the ring 
reinforces the gravitational attraction of the star on the dust layer at the outer edge of the ring, 
and has the opposite effect at the inner edge of the ring. These gravitational perturbations must 
be balanced by pressure gradients or centrifugal forces. Since H p < Ar < H g , the orbital speed 
of the gas is much less affected than the orbital speed of the dust, which increases at the outer 
edge of the ring and decreases at the inner edge. The outer (inner) edge suffers greater (less) drag 
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and a larger (smaller) inward radial drift rate. This reinforces the perturbation in S p . Drag is 
clearly essential: without it, contraction of the dense ring would be prevented for Ar 3> H p by the 
difference in specific angular momentum between the inner and outer edges of the ring; but drag 
changes the angular momentum of the dust. The planform of the secular instability is expected to 
be stationary in a frame drifting inward with the dust layer. Safronov (1991) remarked that the 
layer will probably be stabilized by the difference in radial drift speed among bodies of different 
size. This assumes that the drag acts separately on individual bodies. Many theoretical studies of 
the dust layer (including GW and CDC), however, assume that the drag acts collectively because 
the gas within the layer is entrained by the greater mass of solids. In §4, we shall find that whether 
the particles feel the drag individually or collectively depends upon their size. For the present, we 
take the drag to be collective. 

The original intent of the present work was to study the secular gravitational instability de- 
scribed above. We expected to find a slowly growing mode (compared to the orbital frequency) 
on scales Ar » H p , whose growth rate would tend to zero with p/p*. But our analysis uncovered 
a rapid, small-scale, drag-driven instability that does not require any selfgravity, although it does 
require that the drag be collective. 

The structure of this paper is as follows. In §2, we present height-integrated nonlinear equa- 
tions for the axisymmetric evolution of the dust layer based largely upon the work of CDC, and 
we cast these equations into dimensionless form. We also discuss our assumptions for the nebular 
parameters, and we find steady solutions to the equations in which the surface density is locally 
uniform. The linear stability of the steady solutions is examined in §3, first analytically and then 
semi-numerically. The most rapid instabilities are driven by drag rather than selfgravity. The fol- 
lowing section, §4, explores a number of important issues, including the collective nature of the drag 
force, the relative importance of the instability for building planetesimals compared to straightfor- 
ward collisional agglomeration, and the expected planetesimal mass if the instability dominates. 
§4 also briefly compares this work with some notable previous investigations of the stability of 
protoplanetary disks that treat gas and dust as dynamically distinct components (Coradini et al. 
1981; Sekiya & Nakagawa 1988; Noh et al. 1991). §5 sums up. 



2. PHYSICAL MODEL AND GOVERNING EQUATIONS 

The purpose of this section is to develop model equations for axisymmetric radial perturbations 
of the dust layer. To avoid having to treat two spatial dimensions at once, we adopt a height- 
integrated approximation based upon the results of CDC, whose notation we generally follow. 
Important physical influences on the evolution of the dust layer include turbulent drag against the 
gas; turbulent mixing, horizontally as well as vertically; gas pressure (which resists compression of 
the layer); and selfgravity. 

We are concerned with a late phase in the evolution of the gaseous protoplanetary nebula when 
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accretion and its attendant turbulence have subsided, allowing the dust to settle into a thin layer 
of half-thickness H p very much less than that of the gas, H g (Safronov 1969; Goldreich & Ward 
1973). (In this paper, "dust" refers indiscriminantly to all solids smaller than a few meters.) Within 
the dust layer, the volume-averaged mass density of the dust (p p ) exceeds that of the gas (p g ), as 
shown below, but the column density of the gas (S g « 2H g p g ) greatly exceeds that of the dust 
(S p 2H p pp). At solar abundance, S g /S p « 200 if the dust contains all materials condensible at 
temperatures ~ 300 K. 

We assume a standard minimum-mass solar nebula with the following radial scalings: 

-0.5 



T(r) = 280 K(^)- 



S g (r) = 1700 g cm" 2 ( —\ ^ 



VAU/ 



/ r \ -2.75 

E (r) = 1.7 x 10- 9 g cm-3 (_) . (2) 



For ease of reference, we summarize the more important physical parameters in Table I. Table 

I 

here 

2.1. Turbulent drag on the dust layer 



In steady state, the dust layer orbits nearly at the the keplerian velocity Vk = \J GM/r, where 
M is the mass of the protostar, and r is the orbital radius. (Unless otherwise stated, M = 1 Mq 
and r = 1AU.) Because of its radial pressure gradients, the gas orbits at the slightly lower velocity 
(1 — t])Vk determined by equation (1). For nominal parameters, rj ~ 10~ 3 [Table I]. 

Because of the orbital velocity difference above, GW argued that the interface between the dust 
layer and the relatively pristine gas should be turbulent, even if the nebula is otherwise quiescent, 
by analogy with the standard laboratory situation of a thin solid plate placed edgewise in a nearly 
inviscid, incompressible flow. The drag on such a plate due to turbulent momentum transport is 



Pel Aw I . . 

^drag = ~% l Av, (3) 

-tie* 



where Av is the velocity difference between the plate and the undisturbed flow, and the "Reynolds 
number" Re* is an empirical dimensionless constant. The drag law is usually stated in scalar 
form, .Fdrag = — Pg(Av) 2 /Re*; the vectorial form above assumes that the force is antiparallel to 
the relative velocity. Citing laboratory measurements, GW took Re* ~ 500, but CDC advocated 
Re* f» 45 — 180. We take Re* = 180 as a reference value but are careful to note how things scale 
with this number. 

The turbulent transfer of momentum from the dust layer (or flat plate) to the surrounding 
medium is effectively diffusive and can be described by a turbulent viscosity ut, with dimensions 
of (length) 2 (time) -1 . In order to reproduce the drag (3), one needs is? w LA-u/Re*, where L is 
the vertical thickness of the turbulent boundary layer — the layer within which the velocity, though 
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fluctuating, is closer to that of the dust layer than to that of the undisturbed gas. As noted by 
GW, L should be comparable to the Ekman length constructed upon is? and the orbital angular 
velocity VL: 

Viewed in inertial coordinates, the direction of the stress exerted by the dust layer on the gas 
rotates during the orbit, and Le is the vertical distance over which horizontal momentum diffuses 
in time fi -1 . Similar considerations govern the thickness of the boundary layer at the base of the 
Earth's atmosphere (Pedlosky 1987). As a generic measure of the boundary-layer thickness, CDC 
take L = c^Le with cl = 1.5. Since u T rs LAv/Re*, it follows that 

K -fez- (5) 



and also that 



c L Av 

ReTFT (6) 



The factors involving Re* in eqs. (3), (5), and (6) represent significant corrections to naive 
dimensional analysis. They are semi-empirical, based upon terrestrial experience that may not be 
a reliable guide to the situation we wish to study. However, it seems difficult to improve upon these 
factors without resort to an ambitious three-dimensional hydrodynamic simulation. 



2.2. Thickness of the dust layer 

The equilibrium half-thickness H Pt o of the dust layer results from competition between turbu- 
lent stirring and gravitational settling. Following CDC, we assume that H p is comparable to the 
1% boundary-layer thickness — the height above the midplane at which the mean azimuthal velocity 
is within 1% of its asymptotic value far above the dust. Also following CDC, we take this thickness 
to be larger than L by a factor c$ ~ 2.5. Thus, 

Unless otherwise noted, we adopt eq. (7) for the dust-layer thickness. We believe that it is appro- 
priate for dust-particle sizes of order 10 cm. But the layer is likely to be somewhat thicker (thinner) 
if the particles are much smaller (larger) than this. 

When all particles are small they will be tightly coupled to the gas within the dust layer by drag 
(see below). In this limit, where the settling time is effectively infinite and the slightest turbulence 
is sufficient to loft the dust, the equilibrium vertical distribution of dust is likely to occur at the 
point of marginal shear instability, absent other sources of turbulence. The following is a simplified 
estimate of the thickness under these assumptions; a more careful treatment has been given by 
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Sekiya (1999). Shear instability is expected only if the Richardson number 

R . ^ (d$/dz)(d]np/dz) 1 
(dv/dz) 2 ~ 4" 

The vertical gravitational potential $(z) is Q 2 z 2 /2, and the horizontal velocity v(z) (relative to the 
keplerian value Vk) is determined by the geostrophic balance 

dP 

2pQv = — . (9) 

For simplicity, we take gaussian profiles for both the dust and the gas, p p (z) = p p (0)e~ z ^ 2H p and 
p g (z) = /9 g (0)e~ z2//2H s, with H p < H g . Then since v(z) — »■ -77 Vk when z > H p , 

p(z) = p g { z )+p p { z )^p(Q)(l-f + fe- z2 l 2H l) 

v(z) « + e -*W ~* 

where / = /9 p (0)/p(0) is the fractional contribution of the dust to the total density at the midplane. 
With these assumptions, Ri(z) ~ Ri(0) within one scale height of the midplane, and increases at 
higher altitudes. 

ffi(0) -^M 2 a-/) 2 

On the other hand, 



VV K J P 



Pp (0)H p _ S p 
Pg (0)F g S g p ' 

where Z p w 0.005 is the fractional abundance of dust. Together with H g = c s /£l ps ^^Vk, these 
relations imply that the condition Ri(0) < 1/4 is equivalent to 

p p ( 0)> /4^y/^ 



P (o) ~ V v 

i.e. about equal contributions to the midplane density from dust and gas. It follows that the 
dust-layer thickness in the limit of very small particles is 

(H p )' = Z p H g < 0.18^ (10) 

This is about 9 times larger than our standard estimate (7) if we take Re* = 180 and Av = t]Vk. 
In the small-particle limit where (10) holds, there may be very little drag on the dust layer, since 
the turbulence is only incipient. 

Larger particles are imperfectly coupled to the gas, and significant turbulence will be required 
to keep them from settling. A dimensionless measure of the strength of drag coupling is Qt s , where 
the "stopping time" t s is the decay time of the velocity of an isolated particle relative to the gas 
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if subject to no forces other than drag. As discussed by CDC and others, there are at least three 
distinct drag regimes of interest, distinguished by the relative sizes of the the particle radius (r p ), 
the molecular mean- free path (A), the sound speed (c), and the relative velocity between particle 
and gas (t> re i). Particles with radii r p < 9A/4 are in the "Epstein" regime, where the velocities of 
gas molecules striking the particle are uncorrelated with that of the particle itself, and the drag 
force is 4Trr p p g cv rc \/3. If r p > 9A/4 but r p < Ac/(4f rc i), there is a laminar viscous flow over the 
particle, and the drag is given by Stokes' formula 67r/ir p v Te i, where p « p g \c/2 is the dynamical 
viscosity of the gas. At still larger values of r p and/or v TC \, the flow over the particle is turbulent 
and the drag force increases faster than linearly with v re \, eventually quadratically. The mean- free 
path A ~ l(r/AU) 2 cm at the midplane of the nebula (2). Thus 

t _ f 1.2 X 10 4 (r p /cm)(r/AU) 3 s "Epstein", 

* s ~l 5.3 x 10 3 (r p /cm) 2 (r/AU) 2 - 75 s "Stokes", ^ ' 



with the dividing line between the two regimes at r p ~ 2.3(r/AU) 2 ' 75 cm. Presuming that v TC \ < 7]Av, 
the turbulent regime requires r p > 10(r/AU) 2,5 cm. 

Since the vertical gravitational acceleration is Q 2 z, a particle settling through quiescent gas 
reaches a terminal velocity Q 2 z/t s provided £lt s < 1, so that z decays exponentially at the rate 

7" 1 = (n 2 *.)- 1 . (12) 

Turbulence will resist this settling. Following CDC, we describe turbulent stirring by a mass- 
diffusion coefficient 

- VT /I o\ 

^mass — „ • {^3) 

be 

where the "Schmidt number" Sc > 1 describes the efficiency with which turbulent eddies are able to 
pick up the particles; Sc increases with particle size and with orbital radius in the nebula [cf. Fig. 
2 in CDC]. Therefore, still another estimate of the dust-layer thickness, H£, follows from equating 
the diffusion rate v mass /H 2 to the settling rate Q 2 t s : 

TT,f _( V ^ _ H P (U ) 

p \ScQHJ csVS^m' { ' 

For 10-cm particles at 1 AU, Sc ~ 3 and the thicknesses (7) and (14) are comparable. For much 
smaller particles, where H 1 ' > H p , the turbulence and its associated drag are probably incompletely 
developed, so that the estimate (10) is probably more accurate. For larger particles or larger 
orbital radii, H' p ' is smaller than H p , so that in using eq. (7) for the thickness, we are probably 
underestimating the effects of selfgravity. 



2.3. Vertically integrated equations of motion 



Here we lay out the equations of motion for the dust layer in a height-averaged approximation. 
The gas above and below the dust layer is taken to be undisturbed, on the grounds that we are 
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concerned with radial scales Ar that are somewhat greater than H p but <C H g . Disturbances in 
the pressure and density of overlying gas are expected to equilibrate on a timescale 

c H g 

whereas we shall be concerned with motions of the dust layer over times > f2 _1 . We shall allow, 
however, for the pressure of the gas trapped within the dust layer. 

We write U and V + Vr for the vertically-averaged radial and azimuthal velocities of the dust 
layer: 



oo 

£ p (r,i) = J p p (r,t,z)dz 



V(r, t) J E p (r, t) J rpy ' ' ' \ v^r, t, z) - V K {r 



Continuity of dust mass is expressed by 

as 



d£ p 1 d 
— - H 

dt r dr 



r I UY, p - v—^ 
or 



0. 



The term involving v reflects our assumption that if the turbulence mixes the particles vertically, 
then it also mixes them radially. For want of a better hypothesis, WG tclke V — ^mass 

[eq. (13]. On 

lengthscales Ar <C r, the metric factors r _1 and r surrounding the radial derivative above can be 
neglected, so our mass equation becomes 

Next, we develop equations for the velocity components V and U. In axisymmetry the only 
azimuthal force is drag, so that if radial diffusion is neglected for the moment, then the specific 
angular momentum of the dust evolves as 

1 + ^1)^ + ^)1— r^(V + ,l4), 

or equivalently, since V/r <C VL, 



dt dr J 2 SpRe* 

where 

Av = \Av\ = ^U 2 + (V + 7]V K ) 2 . (17) 



-9- 



The dynamical equation for U incorporates coriolis force (the residual of centrifigal and central 
forces), drag, self-gravity, and pressure: 

An adequate local approximation for the acceleration due to self-gravity is 

oo 

- a >^J^TW^ d{ ^ (19> 

— oo 

which is equivalent to its radial fourier transform, 

oo 

tp(k) = fi/>(r + Ar) e~ ikAr d(Ar) = e -l fc l H p° E(fc), (20) 

J 1^1 

oo 

if the main contribution comes from |Ar| <C r. The terms involving -ff p ,o soften the force on 
scales smaller than the equilibrium dust-layer thickness. Perhaps one ought to use the actual local 
thickness H p (r,t) rather than its equilibrium value, but then the force could not be calculated so 
easily by fourier transforms. It seems pointless to accept this complication since the exact force 
softening depends upon the detailed shape of the vertical density profile, which is beyond the scope 
of the height-integrated treatment. 

The height-integrated pressure involves the weight of the overlying dust: 

oo oo 

II(r,t) = J dz j dz p p (r,t,z f )n 2 z' = £ p Fp^ 2 , ( 2 1) 

OO 2 

neglecting a factor of order unity that also depends upon the vertical density profile. This formula 
presupposes vertical hydrostatic equilibrium, which is justified by the thinness of the dust layer 
and the slowness of its motions. Only the dust density contributes to the excess pressure since the 
gas is neutrally buoyant. 

We assume that the gas within the dust layer is trapped there on timescales shorter than the 
settling time (12). The internal (thermal) energy of this gas is large compared to the kinetic energy 
in the radial motions of interest, even if the trapped gas is a minor contributor to the volume mass 
density. On short timescales, therefore, the dust layer is is three-dimensionally incompressible: 

fp.fp and s3>. (22) 

2j p Hp II Zj p 

But Hp should relax towards equilibrium on the settling time 7" 1 [eq. (12)], restoring II to its 
equilibrium value. These considerations suggest the following dynamical equation for II: 

s + c, l) n = - 3n f -^< n - n «>. < 23 > 
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where IIo is the equilibrium pressure. 

For the purpose of nonlinear simulations, it is useful to re-cast the dynamical equations (16), 
(18), and (22) in "conservative" form. Multiplying eq. (16) by S p and eq. (15) by V and adding 
the results, we have 

at or dr A 2 Re* 

The form of the term involving v results from our assumption that mass diffuses radially but 
velocities (momentum per unit mass) do not. This is unreasonable, since one already assumes 
that vertical momentum diffusion is driven primarily by vertical shear. It is very unclear, however, 
exactly what combination of radial derivatives of V and S p to use above. There is a controversy 
whether hydrodynamic turbulence in a differentially rotating disk attempts to equalize angular 
velocity, (Vk + V)/r, or specific angular momentum, t(Vk + V) (Balbus & Hawley 1998). Without 
taking a stand on this issue, we replace the diffusive term above with 

dr 

for the form of the diffusive term, and corresponding forms for the other dynamical variables. This 
leads to somewhat simpler equations than other plausible choices. More importantly, the radial 
diffusion is then unambiguously stabilizing at short wavelengths; this is a nontrivial result, because 
some other plausible forms of the effective viscosity actually cause axisymmetric instability (Schmit 
& Tscharnuter 1995). 

Our final set of dynamical equations is therefore 



dt dr V dr 



Re* 



(V + r)V K ), (24) 



dt dr V dr 



3X P 
dt 


+ 


d 
dr 


du 


+ 


d 


~dl 


dr 



UZ, - ^) = 0, (26) 

t/n-^)=-2n^-2 7s (n-n ), (27) 



which must be supplemented by Poisson's equation (20). 
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2.4. Dimensionless units 

It is convenient to adopt the system of units and dimensionless variables described in Table II, in 
which ro is a radius characteristic of the region of interest. The numerical values in the final column Table 
of the Table were computed for a minimum- mass solar nebula in "equilibrium" with ro = 1AU, II 
r p = 10cm, Re* = 180, and Sc = 3. In these dimensionless units, equations (24)-(27) and (20) 9 oes 
become 

9(E.V.) + a / w> _^^\ 1 a (28) 



here 



dt* <9x* V cte* 7 2 2 



at* ax* V ax* / 2 ax* ax* 

(29) 

+ flZ" ^-^aZ- = °> (3°) 



:n 1 : + A (V*n* - z,*fM = -2n*f^ - 2 7 *(n* - n* ), (31) 



ar ax* \ ax* / ax 

^(fc*,t*) = 2 —e-\ k *\ H *»t{K,U). (32) 

We have used the abbreviation 

R* = VWTlT+v^ (33) 

for the dimensionless form of the relative velocity (17). 



2.5. Constant states 

We will want to perform a linear stability analysis of eqs. (28)-(32) for small departures from a 
background situation in which the dependent variables are constant with respect to both x* and t*. 
We call this a "constant state." When all of the derivative terms are set to zero, we have II* — > II*o, 
if) if — > an irrelevant constant, and two algebraic relations, 

l^ + (, + v.4 

This leaves a one-parameter family of possible choices for (V*, £/*,£*). Astrophysically, the choice 
is determined by the dust-to-gas ratio; at solar abundance, this is ~ 1/200, leading to £*o ~ 60 at 
1 AU as shown in Table II. While the global ratio is fixed by the initial conditions, the local one 
may vary, since even in the constant state, the dust drifts radially inward with respect to the gas. 
Depending upon the radial profile of gas density and temperature (which we have approximated 
by constants in our local approximation) , this may lead to an enrichment or dilution of the dust at 
small radii. 
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So it is interesting to consider, at least briefly, the full range of possible constant states. The 

simplest equations result if one expresses E* and U* in terms of V*. The physically allowed range 
is 

-K V, < 0, (34) 

and 

17* = -yJ-W^l + V*), 

s * = ( 1 + V *V^T" (35) 

The dimensionless surface density is a monotonic function of V* , but the radial velocity £/* is most 
negative at V* = —1/2, where — ?7* = 1 = £*. (In physical units, this would imply an orbital decay 
time ro/|J7| of only 160 yr at 1 AU!) In the limit V* — > —1 and — ► oo, drag is negligible because 
of the large inertia of the dust, which then follows exact keplerian orbits. In the opposite limit 
V* — ► and — > 0, the inertia of the dust is negligible and drag locks it to the gas. In both cases, 
the radial drift velocity vanishes. The simplicity of equations (35) motivated the choice of units 
shown in Table II, even though the nominal astrophysical values of the dimensionless variables are 
far from unity in this system. 

A very important quantity is the radial mass flux, 

F* = £*£/* = -V(1 + K) 3 (1-3V*) « -1 + 0(V?) if |K| < 1, (36) 

so the mass flux is almost independent of the dust surface density in the astrophysically relevant 
parameter regime. Mass is the only globally conserved quantity in the dust layer — angular momen- 
tum is "lost" to the gas — and we may imagine that F* is fixed by a boundary condition at large 
radius. The instabilities found below may be related to the fact that such a boundary condition 
constrains the constant state only very weakly. 

3. LINEAR STABILITY 

§3.1 analyzes a simplified problem in which both the self-gravity and the pressure of the dust 
layer are ignored. These simplifications allow the dispersion relation to be presented relatively 
easily in closed form. §3.2 presents numerical growth rates that incorporate both of these effects. 

3.1. Analytic treatment without selfgravity or pressure 

The simplified analysis is interesting for several reasons. First, it shows that self- gravity 
is not required for instability, the primary cause of which is the competition between drag and 
inertial forces. Second, the analysis provides a check on our numerically computed growth rates. 
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Finally, there are parameter regimes in which the neglected effects are unimportant — e.g., pressure 
is inconsequential for large particles, which settle quickly to the midplane. 

With II* and Y>* neglected, only three of equations (28)-(32) are needed. We linearize the 
remaining three about an arbitrary constant state (35) and assume that the first-order quantities 
have the dependence 

(Ac/*, A 14, AS*) oc exp [(z — v*k 2 )t* + ik*x*\ . (37) 

After considerable algebra (aided by commercial symbolic-manipulation software) we obtain a cubic 
dispersion relation for the shifted complex growth rate z: 

= [(l + 4w 2 ){l + w 2 ) 3 ] z 3 + [3w(w 2 + l-2ih){l + 4w 2 )(l + w 2 ) 2 ] z 2 
- [(A8h 2 w 4 + \2k 2 w 2 - 8w 8 - I9w 6 - 15u> 4 - 5w 2 - 1 + 56iw 6 h 

+70iw 4 h + Uiw 2 h)(l + w 2 )] z (38) 
+ 8ihw 3 (-6w 4 + 8iw A k if + 2%^ + Ak 2 w 2 - 3w 2 + k 2 - 3u> 6 + Wiw 2 h) . 

To avoid radicals in the coefficients, we have introduced an auxiliary variable w defined by 

w 2 

V* = ? , < w < oo. (39) 

1 + w A 

Nevertheless, the dispersion relation is cumbersome and best understood in limiting cases. 

First, as!i)-> oo, i.e. the state (£/*, 14, £*) — ► (0, —1, 0) where the dust layer is pinned to the 
gas, there are two rapidly damped roots z — w and z ps — 2w. The interesting root is the third 
one: 

z 3ik*w~ l + -klw~ 3 , {w -> oo), (40) 

where we show only the leading order terms in the real and imaginary parts. The actual growth 
rate, however, is Real(z) — v*k 2 . Therefore, since X* pa w~ 2 in this limit, the system is stable if 
S* < (2I/*) 2 / 3 . 

The more important limit is w — > 0, (C/*, 14, £*) — > (0, 0, oo). If one assumes that fc* and z are 
of order unity and sets w = above, then the dispersion relation (38) reduces to z 3 + z = 0. This 
has two purely imaginary roots z = ±i, corresponding to modes that will damp for any v* > 0. To 
study the root near z«0, one has to go to higher order: 

z w (2k 2 - ik 3 )^- 3 [w 0, ->■ oo, fc* ~ 0(1)]. (41) 

This indicates instability if £* < (2/iz*) 1 / 3 . In extracting the approximate root (41) from eq. (38), 
we assumed that fc* is of order unity. But if we take &:* = q/w with g ~ O(l) as w — >• 0, a different 
ordering of terms results, and the dispersion relation reduces approximately to 

(to') 3 - (J) + 2q = 0, u/ = iz + 2g. (42) 
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Clearly, if \q\ is sufficiently large, then there are a pair of complex-conjugate roots for u', implying 
the existence of a root for z with positive real part. The corresponding growth rate is approximately 

31/22-2/3^/3 _ ^(^2. 

The maximum growth rate is therefore 

Real(z) max « 0.49£; 2/ V V5 (£* 00). (43) 
which is achieved at q « 0.31(S^*) -3 / 5 . But (42) has complex roots only if |g| > 3 -3 / 2 , and so 

— 1/2 

this growing mode exists only if £* < 1.5i/# , approximately. 

Scrutiny of Table II shows that the combination S*^* is independent of the critical Reynolds 
number Re*. Therefore, the maximum growth rate is formally independent of this most uncertain 
parameter. The same is true for the ratio 27r(/c*.ff*o) _1 of wavelength to dust-layer thickness at 
maximum growth. 

2/3 1/2 

To summarize, if < <C 1, then constant states are unstable when vj < £* < 
Selfgravity will probably enlarge the range of instability, while pressure will probably reduce it, but 
both effects introduce additional uncertain parameters. 



3.2. Numerical Growth Rates 

As a complement to the analytic treatment of the previous section, we performed a numerical 
linear stability analysis of the governing equations so as to investigate the effects of pressure and 
self-gravity upon the expected growth rates. We began our analysis by transforming the governing 
equations into Fourier space. In this way, the gravitational potential can be eliminated in favor 
of the surface density via Poisson's eq. (20). We next recast the equations by introducing new 
variables for the radial and azimuthal mass flux; / = £*?7* : j = £*V*. We then proceeded 
to linearize the governing equations about the constant states by adding a perturbation to each 
dynamic variable (ie. / — > / + A/), expanding, and then only retaining those terms which arc 
linear in the perturbations. Thereby, with z = d/dt* — ^*<9 2 /<9x*, we wrote the system of linearized 
governing equations in the form 

z(y) = Cy, 

where y = (A/, Aj, AS*, AIT*) is the vector of first-order dynamic variables. Then, given the form 
of the (x*,t*) dependence (37), it follows that the eigenvalues of the linear operator matrix, \c, 
correspond to dimensionless growth rates = Re(A,c) - v^k 2 . 

This analysis was performed for a range of plausible nebula parameters, specifically; r (or- 
bital radius) G {0.1 AU, 1.0 AU, 10 AU}, r p (particle size) G {0 cm, 10 cm, 100 cm}, Re* G 
{45,70,180,500}. Figure 1 shows the growth rate as a function of wavenumber for various pa- 
rameter values. For reference, the unperturbed surface density in our model is about an order 
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of magnitude lower than the critical surface density required for the Goldreich-Ward dynamical 
instability. Fig. 1 

For most choices of the background parameters, the growth rate is largest at the shortest 
available wavelengths. We present results down to a wavelength comparable to the dust layer 
thickness, as the validity of the height-integrated approach is questionable beyond this point. There 
are a few sets of nebula parameters, however, for which the peak growth rate occurs on a significantly 
larger spatial scale. Fig. 2 shows the growth rates for nebula parameters r = 0.1 AU, r p = 10 cm, 
Re* = 180. Again, rates are shown with and without self-gravity. These growth rates are notably 
different from those at other parameter values in that they are considerably lower and become 
negative for all wavenumbers in the absence of self-gravity. 

It is our interpretation that there are two distinct secular instabilities in this system. One is 
the self-gravitating secular instability discussed in §1. It has a relatively low growth rate and a 
long radial wavelength compared to the dust-layer thickness. This is the instability which appears 
in Fig. 2. The other is the drag-induced instability discussed in the previous section, which has Fig. 
larger growth rates and shorter wavelengths, and is exhibited in Fig. 1. Probably because of its 2 
large growth rate, the drag instability is found over a wider range of relevant nebular and particle 
parameters. The rest of the results in this section refer to instances in which the drag instability 
dominates. 

Table III. presents the peak growth rate achieved for a number of other prominent sets of 
nebula parameters. Generally, the dimensionless growth rates Table 

III 

i) increase monotonically with increasing orbital radius. An increase in orbital radius implies ^ ere 
a decrease in the surface dust density and an increase in the dust scale height. [The orbital radius 

also effects the stopping time by determining the drag regime.] 

ii) increase with increasing particle size. An increase in particle size implies an increase in the 
settling rate. 

iii) increase with increasing critical Reynolds number. An increase in Reynolds number implies 
an increase in the dust surface density and a decrease in the dust scale height. 

We repeated the above analysis with the terms relevant to pressure and self-gravity selectively 
neglected. 1 As expected, the growth rates are larger without pressure and smaller without self- 
gravity. But the changes are typically on the order of ten to twenty percent. When both pressure 
and self-gravity are neglected, we recovere the growth rates predicted by the cubic dispersion 
relation (38). 

In summary, our linear analysis reveals a drag-induced instability which grows most rapidly 



1 Recall that in this context, "pressure" means the excess height-integrated pressure due to the weight of the local 
dust column, cf. §2.3. 
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on short wavelengths at large orbital radii, for large particles, and high critical Reynolds number. 



We have shown that for plausible nebular parameters, the growth rate of the drag instability 
can be comparable to the orbital frequency fi. 

The instability is rather robust. The following toy model shows that the tendency to instability — 
as opposed to the actual growth rate — does not depend upon the more intricate and uncertain as- 
sumptions that we have made in formulating our dynamical equations. Consider a one-dimensional 
distribution of mass a(x, t) per unit length that is subject to a fixed "gravitational" acceleration 
g and an opposing "frictional" force per unit length —f(a)v. Here ag is a proxy for the gravi- 
tational and coriolis terms in the dust-layer equation of motion — forces that are proportional to 
mass — whereas / stands for the drag force on the surface of the dust layer, which is essentially 
independent of the surface mass density. The equations for mass and velocity (v) are 



This system has equilibria with uniform velocity vq and surface density o"o satisfying o^g = /(co). 
Without loss of generality, vq = 0. Linear perturbations with (x,t) dependence exp(st + ikx) have 
the dispersion relation 



Clearly, there will always be an unstable root, Real(s) > 0, regardless of the details of the friction 
function /(<r), unless the frictional force is strictly proportional to mass. One could certainly 
complicate to this model to stabilize it, but the point is that a balance between inertial and non- 
inertial forces tends to be unstable. 

The one assumption that is truly critical for instability is that drag should act coherently upon 
the dust layer as whole. If the drag acts upon each particle independently, collective instability 
cannot result from drag alone. The assumption of collective drag, introduced by GW, is reasonable 
if the gas within the dust layer is substantially entrained by the dust, which requires that the 
turbulent wakes of particles overlap before they mix with the gas above and below the dust layer. 
The net drag will then be less than it would be if each particle independently encountered the full 
headwind from the subkeplerian gas. Thus, a necessary condition for the drag to be treated as a 
collective effect is that the force per unit area (3) be less than N x f(r p ), where f(r p ) is the drag 
on an isolated (spherical) particle of radius r p , and N is the number of such particles per unit 
area. For simplicity, let the dust layer be composed entirely of particles having this radius; then 
the surface density of the layer is related to r p by 



4. 



DISCUSSION 



d t a + d x (av) = 0, 

d t v + vd x v = g - vf(a)/a, 





(44) 
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where p s « 1 g cm 3 is the density of the solid material. In the size regime r p ~ 10 — 100 cm, the 
drag on this particle can be estimated from Stokes formula, 

/(r p ) = 67rp g r p Av, (45) 

with dynamical viscosity p g « 0.9 x 10~ 4 g/cm-s at 300 K. Thus 

Nf 9Re*S p/ u g ' \180 J VlAU/ V 100 cm/ ' k j 

using the nebular parameters and scalings of equations (2) and Table I, and assuming that AV ~ 
c 2 /Vk- Actually, the force ratio (46) is independent of the temperature and density profile of nebula 
as long as the relative abundance S g /S p of gas and dust is constant. We conclude from this that 
a collective treatment of the drag is appropriate as long as r p < 3 m. 

Several groups have previously analyzed the stability of disks when gas and dust are treated 
as dynamically distinct components, but none of them seem to have identified the drag instability 
that we have emphasized, and it is interesting to ask why not. Coradini et al. (1981) performed 
a stability analysis local to the interior of the dust layer, where the dust was treated as uniformly 
mixed into the gas but enhanced in abundance by sedimentation. Their vertically local treatement 
overlooked the exchange of angular momentum between the layer as a whole and the relatively dust- 
free overlying gas. The approach of Sekiya & Nakagawa (1988) was similar, except that they allowed 
for more general perturbed velocities within the layer (not necessarily parallel to the midplane). 
Both of these groups found instabilities, but only for dust densities larger than the critical value p* 
(cf §1), so these presumably were modes that rely upon selfgravity. Noh et al. (1991) investigated 
global nonaxisymmetric stability of the disk. They allowed for a thin dust layer but represented 
the drag force by a term strictly proportional to mass, which suppresses the nonselfgravitating 
drag instability, as the toy model above demonstrates. Noh et al. (1991)'s approximations should 
have captured the secular selfgravitating mode sketched in § 1 if they had considered axisymmetric 
modes; they may have found a nonaxisymmetric version of the instability but did not interpret 
their results along these lines [see their Fig. (4b)]. 

Is the drag instability necessary to form planetesimals, or will collisional agglomeration may 
accomplish the task more quickly? There are at least two issues: 

• (i) Do collisions result in larger particles, i.e., do colliding particles stick or fragment? 

• (ii) How does the collision rate compare with the growth rate of the drag instability? 

The first question is difficult because one does not understand the structure and elastic strength 
of the particles. If these are loose aggregates, fragmentation is likely at relative velocities > 10 m s _1 
(Weidenschilling & Cuzzi 1993, and references therein). 

As for the second question, it is clear that the collision rate is as fast or probably faster than 
the drag instabilities. If all particles are spheres of radius r p , then the probability of crossing the 
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dust layer without a collision is ~ exp(— 3£p/p p r p ), so that the collision rate per particle 

(«„*)-> ~ ^_ w 0.13 (J*-)" 1 (A)" 3 ' 2 yr - (47 ) 
2p p r p \ 1 m/ V AU/ 

for the nebular parameters of Table I. 

Therefore, if most collisions result in sticking, then collisional agglomeration is probably more 
important than drag instability. But, it is interesting that the latter process is at competitive 
with collisions. The collective instability has the attractive feature that it may take planetesimals 
directly from sizes ~ 10 cm to ~ 1 km in a single step, and thereby bypass the size regime in which 
orbital decay (due to aerodynamic drag on individual planetesimals) is dangerously rapid. 

Assuming that the collective instability does operate, it will be important to understand the 
nonlinear outcome. We have performed time-dependent numerical integrations of the nonlinear 
system (28)- (31), neglecting selfgravity and pressure. Typically, the surface density peaks an order 
of magnitude or more above the background. No simple pattern emerges; multiple peaks develop at 
slightly different radial drift speeds. Numerical stability does not seem to be an issue as long as a 
suitably short time step is used, but unfortunately, the detailed nonlinear outcome depends upon the 
form chosen for the diffusive terms. In eqs. (28)-(31), the diffusive terms are in "flux-conservation" 
form, but there is no clear reason to demand such a form except in the mass equation (31), since 
radial and azimuthal momentum can be exchanged with the gas. We found nonconservative forms 
that produced the same linear growth rates but more organized nonlinear evolution. Because of 
these uncertainties, and because selfgravity and nonaxisymmetry would have to be included in 
order to understand planetesimal formation, we shall not give details of our nonlinear experiments. 

Nevertheless, we cannot resist using the linear results to speculate about the mass of the 
planetesimals that may form by this instability. Since the linear growth rate appears to peak at a 
wavelength of order 4irH p (Fig. 1), a plausible estimate for the mass is 

M p i anetesimal ~ £ p (2^ p ) 2 ~3x 10 19 g. (48) 

The numerical value is evaluated for the surface density and thickness of the dust layer at 1 AU 
(Table I). After expulsion of gas and collapse to solid densities, this corresponds to a planetesimal 
radius ~ 10 km — larger than for GW's direct gravitational instability because our dust layer is 
thicker than theirs. 

5. SUMMARY 

We have found a new collective instability of the dust layer in a quiescent protostellar disk. 
Like the gravitational instability contemplated by Safronov (1969) and Goldreich & Ward (1973), 
the new one may collect particles smaller than a meter directly into planetesimals larger than a 
kilometer, thereby circumventing a possible difficulty with rapid orbital decay of the particles. But 
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the new instability is driven by drag rather than selfgravity. In fact, it is driven by the very shear 
turbulence that limits the density of the dust layer and forestalls direct gravitational instability. We 
envison a two-stage mechanism in which drag first concentrates the dust radially into rings until 
they become selfgravitating, whereupon the rings fragment nonaxisymmetrically. We find that 
instability is generic when inertial forces balance frictional forces in a nonrigid mass layer. Thus, 
the tendency to instability is independent of the details of the dynamical equations we have used 
to describe the dust layer. To the extent that these equations are correct and that the parameters 
of a standard minimum-mass solar nebula are applicable, the growth time of the drag instability is 
of order the orbital period and comparable to the two-body collision time. Planetesimals formed 
by this instability are estimated to have radii of order 10 km at 1 AU. 

These conclusions are subject to several caveats. Most importantly, we lack a fundamental un- 
derstanding of the exchanges of momentum and mass between the dust and gas layers. Turbulence 
is central to these exchanges, and we are unaware of any directly comparable terrestrial analogs 
to the dust layer, which is probably better described as a slurry than the flat plate envisaged 
by Goldreich and Ward. Formally at least, our growth rates are rather insensitive to the critical 
Reynolds number Re*, which describes both the drag on the dust layer and its thickness (§2 & §4), 
and to the particle size. Most likely, however, these parameters must lie within a limited range if 
our dynamical equations are to be valid. It is vital for the drag instability that the drag on the 
dust layer be treated as a collective process. According to our best estimates, a collective model 
is appropriate for particles smaller than about 3 meters at 1 AU, and this limiting size depends 
upon Re* (§4). On the other hand, when the particles are much smaller than 10 cm, shear-driven 
turbulence may not be fully developed, so that the drag and the growth rate may be substantially 
reduced (§2.2). Also, the equilibration time of the dust-layer thickness is a strong function of parti- 
cle size; there may be parameter regimes in which drag quickly produces large contrasts in surface 
density, but gravitational instability has to wait for particles to settle and the volume density to 
increase. Direct numerical simulations of dust-gas mixtures in three dimensions may be necessary 
to answer these questions. 

This work was supported by the NASA Origins program under grant NAG5-8385. 
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TABLE I 

Physical parameters of the protoplanetary disk and dust layer: symbols, meanings, 
and assumed numerical values. 



Symbol 


Meaning 


Formula 


Nominal Value (1AU) 


c 


isothermal sound speed 


y/kuT/ji 


1.0km s" 1 


CL,C«5 


boundary-layer factors 




1.5, 2.5 


V 


azimuthal drift factor 


{V K - V gas )/V K 


10~ 3 


H g 


gas-layer ^-thickness 


c /n 


5.1 x 10 6 km 


Hp 


dust-layer ^-thickness 


CLCsr]r /Re* 


3.1 x 10 3 km 


V K 


keplerian velocity 


VGM/r 


30 km s _1 


n 


orbital angular velocity 




2vryr- 1 


r 


orbital radius 




1 AU 


Re* 


critical Reynolds number 




180 


Sc 


Schmidt number 




1.0 


S g 


gas surface density 




1700 g cm" 2 


z P 


dust abundance by mass 




0.005 




dust surface density 




8.5 g cm -2 


U 


radial drift speed 


-2p g {r]V K ) 2 /Re*nS p 


— 1.0m s _1 
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TABLE II 
Dimensionless variables and parameters. 



Symbol Meaning Definition Nominal Value 



Dimensional Units 

Ti time rr 1 1 yr/2vr 

L\ length rjro 10 -3 AU 

S\ surface density 2p g rjro/Re* 0.28 g cm~ 2 

Independent variables 

t* time t ■ f2 

radial separation (r — ro)/Li 

k* radial wavenumber k • L\ 

Dependent variables 

X* surface density S p /5i 30. 

XJ* radial drift speed U/rjVK —0.033 

V* azimuthal drift V/rjV K -2.8 x 10~ 4 

II* 2D pressure U/S^-nVx) 2 0.013 

ip* self-potential ^/(^^k) 2 

Dimensionless constant parameters 

diffusivity c|/Re 2 Sc 2.3 x 10~ 5 

H^q equil. ^-thickness c^c^/Re* 0.021 

IT*o equil. pressure H*qH% 0.013 

q* self-gravity 7]Vk^I/^GSi 1.0 x 10 4 

7* settling rate Qt s 0.11 
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TABLE III 

Local axisymmetric instabilities of the dust layer. Numerical growth rates (column 
4) are given in units of local orbital frequency Q. 



r p (cm) r (AU) Re* Maximum e— folding Comments 

growth rate time (yr) 






0.1 


180 


1.8 xl0~ 4 


28.0 




10 


0.1 


45 


2.1 xl0~ 4 


24.0 


long-wavelength instability 


10 


0.1 


180 


1.9 xl0~ 4 


26.5 




i n 


n 1 

U.l 


ouu 


U.oU 


1 7 v i n— 2 

1. / X 1U 




100 


0.1 


180 


0.24 


2.1 xl0~ 2 







1.0 


180 


0.46 


0.35 




10 


1.0 


45 


0.45 


0.35 




10 


1.0 


180 


0.55 


0.29 


nominal parameters 


10 


1.0 


500 


0.69 


0.23 




100 


1.0 


180 


0.73 


0.22 







10.0 


180 


0.78 


6.45 




10 


10.0 


180 


0.76 


6.49 




10 


10.0 


180 


0.90 


5.58 




10 


10.0 


180 


1.01 


4.96 




100 


10.0 


180 


0.94 


5.31 
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0.1 0.2 0.3 0.4 0.5 

wavenumber 



FIG. 1. Dimensionless linear growth rates for nebula parameters r = 1 AU, Re* = 180. The four 
curves represent models for r p = 10 cm with (solid line) and without (short dashed) self-gravity, 
as well as for r p = cm(dashed) and r p = 100 cm (dot-dashed), both with self-gravity. The units 
of wavenumber are such that a wavenumber of 1.0 corresponds to an instability at a wavelength of 
2ttH p . 



-25- 



0.001 



0.0005 







o 

&0 



-0.0005 



-0.001 







0.02 



0.04 



0.06 



wavenumber 



FIG. 2. Dimensionless linear growth rates for nebula parameters r = 0.1 AU, r p = 10 cm, Re* = 
180. The two curves represent models with (solid line) and without (short dashed) self-gravity. 
The instability for this set of nebula parameters is slowly-growing, occurs at long wavelengths, and 
is not seen in the absence of self-gravity. The units of wavenumber are such that a wavenumber of 
1.0 corresponds to an instability at a wavelength of 



